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Magnetic three states of matter: A quantum Monte Carlo study of spin liquids 
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We present thermodynamie phase diagrams showing magnetie analog of “three states of matter,” namely, spin 
liquid, paramagnetic, and magnetically ordered phases, obtained by unbiased quantum Monte Carlo simulations. 
Our simulations are carried out for Kitaev’s toric codes in two and three dimensions, i.e., the simplest realizations 
of gapped topological Z 2 spin liquids, extended by a nearest-neighbor ferromagnetic Ising coupling. We find 
that the ordered phase borders on the spin liquid by a discontinuous transition line in three dimensions, while 
it grows continuously from the quantum critical point in two dimensions. In both cases, our results elucidate 
peculiar proximity effects to the nearby spin liquids in the high-temperature paramagnetic phase, even when 
the ground state is magnetically ordered. The thorough study of magnetic three states of matter is achieved by 
introducing the “fictitious vertex” method into the directed loop algorithm. This provides a generic prescription 
to simulate models with off-diagonal multispin interactions, in which the conventional scheme may suffer from 
intrinsic ergodicity breakdown as in the present case. 

PACS numbers: 75.10.Jm, 03.65.Vf, 05.50.+q 

Introduction. The concept of topological order becomes 
central for exploring new states of matter^. Particularly, the 
possibility of quantum spin liquids (QSLs) with topological 
orders has been extensively discussed to explain experiments 
in highly frustrated quantum magnets in which conventional 
symmetry-breaking order seems prevented much below the 
Curie-Weiss temperature^. Naturally, more theoretical inputs 
for verifying QSLs are desired. Whether or not QSLs have 
unmistakable hallmark is one of the most relevant questions. 

Also, even if a system under investigation may develop a mag¬ 
netic order at the lowest temperature, we may ask if there are 
any characteristic proximity effects to a nearby possible QSL 
phase. To answer these questions, the Kitaev modeP’^ and its 
descendants are useful as their ground states are QSLs, pro¬ 
viding good starting points for exploring quantum and thermal 
effects^“^^. Some of them may have connections with real ma¬ 
terials such as iridates after including perturbations stabilizing 
magnetic orders Recently, two of the authors and their 
collaborator investigated thermal effects in the Kitaev models 
on honeycombhyperhoneycomband decorated hon¬ 
eycomb^^ lattices, clarifying the nature of QSL-paramagnetic 
transitions, which is very sensitive to types of QSLs and di¬ 
mensionality. 

A combination of thermal and quantum effects on Kitaev- 
type models is a more challenging topic especially in three 
dimensions (3D). To the best of our knowledge, effects of a 
magnetic field were investigated at both zero and finite tem¬ 
perature (T) in two dimensions (2D)^"^, but not in 3D. As 
for exchange-iypQ perturbations, although they pose interest¬ 
ing questions regarding transitions between topological and 
conventional ordered phases, the corresponding study at finite 
T was not initiated until a recent functional renormalization- 
group calculation in 2D^^. The studies in 3D were so far re¬ 
stricted to semiclassical treatments^^’^^ or the Bethe lattice^"^. 

In this Rapid Communication, we wish to initiate a full 


quantum treatment of the exchange-type interaction effects 
on a 3D QSL at finite T. We also investigate a 2D model 
and compare the phase diagrams. We will focus on gapped 
Z 2 QSLs, which may be realized in certain frustrated Heisen¬ 
berg systems (e.g., 2D kagome^"^ and 3D hyperkagome^^) or 
iridates (e.g., related compounds with Li 2 lr 03 ^^’^^’^^). Since 
our interest is in generic properties of QSLs, we will consider 
Kitaev’s toric codes on square^ and cubic^^“^^ lattices as the 
simplest Hamiltonians realizing the gapped Z 2 QSLs. We add 
a nearest-neighbor ferromagnetic (FM) Ising coupling, which 
introduces quantum fluctuations to otherwise static flux exci¬ 
tations of the QSL states. These quasiparticles are known to 
have nontrivial properties, such as anyonic mutual statistics 
in 2D (Ref. 3) or a conflnement-deconflnement transition in 

3J)19,20,29,30 

Our phase diagram obtained by unbiased quantum Monte 
Carlo (QMC) methods completes a study of the magnetic ana¬ 
log of “three states of matter,” namely, QSL (“liquid”), para¬ 
magnetic (“gas”), and FM (“solid”) phases, while the first two 
were the subjects of Refs. 18-20. We And a qualitative differ¬ 
ence between 2D and 3D and elucidate the distinct nature of 
quasiparticles responsible for the difference. We also unveil 
proximity effects to the QSL phases even when the ground 
state is magnetically ordered. Although conventional QMC 
algorithms turn out to be inefficient in and near the QSL phase, 
or even nonergodic with the directed loop updates^ we 
solve this issue by introducing the “fictitious vertex” method, 
a generic idea to handle off-diagonal multispin interactions in 
QMC in a simple manner. 

Model. Both in 2D and 3D, the Hamiltonian can be for¬ 
mally written as 

s p (ij) 

where S = 1/2 spins [cri = {af ,a/ ,af) are Pauli matri- 
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FIG. 1. (Color online) (a) 3D and (b) 2D lattices for H in Eq. (1) with 
circles representing S' = 1/2 spins, (c) Electric charge As = —1 
(crosses) and magnetic flux Bp = —1 (plaquettes) excitations in 
2D at the ends of strings of and operators, respectively, (d) 
Electric charge excitations in 3D. (e) and (f) Examples of magnetic 
flux loops in 3D with crosses indicating plaquettes with Bp = —1. 

ces] reside on links and (ij) denotes nearest neighbors. A^, 
Xb, and Jxx are positive. As = Ylj^s^j denotes a four- 
(six-)spin product of assigned to a node s of the square 
(cubic) lattice, while Bp = Ylj^p cr^ is a four-spin plaque- 
tte operator [Figs. 1(a) and 1(b)]. Since [As, Bp] = 0 Vs 
and Vp, both in 2D and 3D, these are local conserved quan¬ 
tities of ^ = — Aa ^s Xfs — Xb ^pBp with eigenvalues 
± 1 , though not all of them are independent^’^^. is ex¬ 
actly solvable and the condition for a ground state is that all 
of these eigenvalues are +1. In addition, there are nonlocal 
Z 2 operators [global string operators in 2D (Ref. 3) and ei¬ 
ther global string or global surface operators in 3D (Ref. 29)] 
commuting with any ^4^ and Bp, and thus with They 
distinguish 2 ^-fold degeneracy of the spectrum of on the 
(i-dimensional torus^’^^. The FM Ising coupling Jxx makes 
Bp's unconserved (see below). For small Jxx^ splitting of the 
topological ground-state degeneracy is exponentially small in 
the system size L (Ref. 3). 

In both dimensions, “electric charges,” or nodes with As = 
— 1, can be created or annihilated in pairs by erf (i G s) and 
they are deconfined quasiparticles [Figs. 1(c) and 1(d)]. A 
crucial difference between 2D and 3D appears in “magnetic 
fluxes” corresponding to plaquettes with Bp = — 1 , created 
or annihilated by af with i £ p [Figs. 1(c), 1(e), and 1(f)]. 
While fluxes in 2D are pointlike and deconflned, they always 
form loops in 3D as the operator identity Ilpecube Bp = 1 
implies that a flux string entering an elementary cube must 
leave it (no magnetic monopole). To stretch a flux loop, ad¬ 
ditional plaquettes must be flipped with an extra energy linear 
in the increased length [Fig. 1(e)]. Consequently, loops are 
conflned at low T. More delocalized loops appear at higher 
T favored by entropic effects, suggesting that a proliferation 
of flux loops may occur at a flnite T. In fact, it was shown 



FIG. 2. (Color online) (a) Phase diagram in 3D for A a = Xb- Some 
error bars are smaller than the symbol size. The lines are guides to 
the eye. The intensity plot, presented for T/Xb > 1, corresponds 
to {w%) fovL — A [see also Fig. 3(c)]. (b) Phase diagram in 2D 
for A A = Xb. The dashed line shows the asymptotic behavior of 
T* = 0.2 A ~ O.SXb — 1.6 Jxx, where the prefactor (= 0.2) is 
chosen to match T* with the peak position of Xxx- 


for Jxx = 0 that this corresponds to a flnite-T transition^^’^^. 
Topological entanglement entropy is flnite for T < al¬ 
though expectation values of the nonlocal Z 2 topological or¬ 
der parameters vanish for any flnite 

This transition is expected to be of the inverted 3D Ising 
universality class. This can be explained by the duality 
between the 3D classical Ising model on the simple cubic 
lattice and the flux sector of the 3D toric code. By us¬ 
ing the high-T expansion, the partition function of the Ising 
model can be rewritten as a sum over loop conflgurations: 
Rising oc where K and ^otai are the di¬ 

mensionless coupling and the total loop length of a conflgura- 
tion, respectively. Thus, these models can be related through 
exp(—2A B IT) = tanh K, with a minor caveat that the wind¬ 
ing number of flux loops (high-T graphs) is even (can be 
even or odd), which however would not alter the universality 
classConsequently, the two models should undergo a phase 
transition of the same (i.e., 3D Ising) universality class. How¬ 
ever, this is with the inverted T axes in a sense that conflned 
and extended loops are favored, respectively, in the low- and 
high-T (high- and low-T) phases in the 3D toric code (Ising 
model). Hence, universal amplitude ratios of one model are 
inverses of the other. 

Phase diagram in 3D. For Jxx 7 ^ 0, the flux loops acquire 
quantum dynamics, while electric charges remain static. The 
model is no longer exactly solvable, and we perform continu- 
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FIG. 3. (Color online) FSS plots of (a) {Cb) for Jxxj^B = 0.1 and 

(b) Xxx for Jxxj^B — 0.35, respectively, produced by assuming 
Df = 1.7349(65)^^ ry = 0.036 39(15), and v = 0.630 12(16)^^ 

(c) {vj\) for L = 4. (d) The histogram of the energy density for 
Jxx/^B = 0.3 slIT/Xb = 1.7. (e) for JxxIXb = 0.2, 0.3, 0.4, 
and 0.5. All the data correspond to the 3D model with A a = Xb- 


ous time world-line QMC simulations. We show the obtained 
phase diagram in 3D in Fig. 2(a), leaving the details of our 
method for later. Below, we will discuss the case of Aa = A^, 
where fluxes and electric charges have comparable excitation 
energies for small Jxx • 

The paramagnetic-QSL transition line cannot be charac¬ 
terized by any local order parameter. Instead, we study the 
flux loop proliferation by evaluating the expectation value of 
length {jCb) of the largest loop. We confirm the inverted 
3D Ising universality by finite-size scaling (FSS) analysis of 
(Cb), which leads to excellent data collapse [Fig. 3(a)]. Here, 
we estimate Tc by using the Bayesian FSS method^'^’^^ as¬ 
suming u = 0.630 12(16)^^ and the fractal dimension Dp = 
1.7349(65) estimated for the high-T graphs in the 3D Ising 
modeP^. We And that Tc decreases only slightly with Jxx, 
suggesting that the effective flux-loop tension is weakly renor¬ 
malized by Jxx- As a different measure of delocalization 
of flux loops, we calculate {w^) defined by the summation 
over the squared winding number of each flux loop C, i.e., 
= Ec §c '^here denotes 

the unit vector along the /i axis, is expected to be zero 
(nonzero) for T <Tc(T > Tc) in the thermodynamic limit^^. 
In the QSL regime (Jxx/^b ^ 0.24), we And {w‘^) decreas¬ 
ing continuously and monotonically as lowering T [Fig. 3(c)], 
suggesting that the paramagnetic-QSL transition is always of 
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FIG. 4. (Color online) (a) Estimated ground-state energy in 3D, 
compared with the results of the second-order perturbative expan¬ 
sions (PEs). (b) ESS plot of the energy gap A and (c) Xxx to an¬ 
alyze quantum critical behavior in 2D assuming 77 = 0.03639(15), 
u = 0.63012(16)^^ and z = 1. 


second order in this model. 

The FM phase with (cr^) 7 ^ 0 is stabilized for large Jxx- 
For sufficiently large Jxx, the flnite-T transition is in the 3D 
Ising universality class as confirmed by FSS analysis of the x- 
component magnetic susceptibility Xxx [Fig- 3(b)]. The tran¬ 
sition changes from second to first order for small Jxx, as indi¬ 
cated by the double peaks in the energy histogram in Fig. 3(d). 
This suggests that there is a tricritical point on 

the paramagnetic-FM phase boundary [Fig. 2(a)]. 

On the FM side of the phase diagram, shows a dip 
enhanced near the QSL-FM phase boundary [Fig. 3(c)]. Here, 
{w ‘^) has a sizable magnitude inside the FM phase as in the 
paramagnetic phase, simply because (cr^) 7 ^ 0 implies a ran¬ 
dom configuration of af. For Jxx > J^x^, the upturn of the 
broad dip is associated with the onset of the critical phenom¬ 
ena of the FM transition. For Jxx < Jxx^, the upturn of the 
pronounced dip is suggested to become a jump in the thermo¬ 
dynamic limit. By comparing with the behavior on the QSL 
side, we conclude that the salient dip in {Wb), implying flux 
loops shrinking with decreasing T above Tc, is a proximity 
effect to the QSL phase. The proximity effect to the QSL 
phase can also be observed in the ^-component magnetic sus¬ 
ceptibility As shown in Fig. 3(e), X;^;^ develops a more 
pronounced peak above T = Tc with decreasing Jxx on the 
FM side. 

The QSL-FM quantum phase transition is strongly first or¬ 
der in this model. Figure 4(a) shows the evolution of en¬ 
ergy density at T/Xb = 4/T (T ^ 0 as T ^ 00 ) from 
both sides of the phase diagram, revealing a level cross¬ 
ing at J^^^^^Xb = 0.232(2). By comparing the data with 
second-order perturbation theory, Tsmaii-= — Aa/3 — A^ — 
J|^/{3Ab) andE’iarge-j,, = -XA/3-4:Jxx->^%/{i8Jxx)for 
small and large Jxx, respectively, we And a good agreement 
indicating that both phases are rather stable up to Jxx = Jxx^^- 

Phase diagram in 2D. The fluxes and electric charges in 
2D are both pointlike (no loop structure) and deconfined^ 
[Fig. 1(c)]. Consequently, infinitesimal thermal fluctuations 
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can destroy the QSL state. At T = 0, the QSL is stable against 
local perturbations because these excitations are gapped. We 
evaluate the gap A by applying the second-moment method^^ 
to (^rC^k=o(^)^k=o(^)) corresponds to a propagator 

of a flux pair for small Jxx- This quantity should be as¬ 
sociated with the lowest-energy gap also near the QSL-FM 
quantum transition, as (cr^) becomes nonzero at this point. 
Our FSS analysis of A suggests a quantum critical point 
(QCP) in the (2-^1 )D Ising universality class at /\b = 
0.1642(3) [Fig. 4(b)]. Near the QCP, our results strongly sug¬ 
gest {As) = 1 as T —> 0 and L ^ oo, providing numerical 
evidence that an effective description of the QCP is that of the 
(2 -f1)D Z 2 gauge theoryconsistent with the gapped spec¬ 
trum in the electric charge sector. This observation motivates 
us to introduce dual Pauli matrices^^ relating Bp = and 
(hp' denotes a site between neighboring pla- 
quettes p, p'), and we obtain two copies of the transverse-fleld 
Ising models. This duality argument supports our observation 
of the (2-1-1)D Ising universality. 

The mapping implies that the primary order parameter at 
the QCP is not (erf) but fractionalized component (pf) 
(while the order parameter at flnite-T transition is simply 
(erf)). Since the two operators comprising belong to the 
decoupled copies in the dual description, the scaling dimen¬ 
sion of (erf) is doubled compared to that of the primary held; 
thus we conclude that the magnetic susceptibility does not fol¬ 
low the standard scaling Xxx ^ but Xxx ^ as 

we confirm in Fig. 4(c). Such unconventional scaling can be 
seen as another nontrivial proximity effect to QSL in 2D, and 
should be observable as Curie-Weiss-like scaling in the quan¬ 
tum critical regime at finite T: Xxx ^ Here we note 

that other quantities, such as the specific heat, exhibit clear 
deviations from the paramagnetic behavior. We also evaluate 
Xxx on the QSL side, where we And a non-divergent peak at 
T ^ 0.2A. This defines crossover temperature T* [Fig. 2(b)], 
below which thermodynamic properties are very similar to the 
ground-state properties. The FM phase is stable at T > 0, as 
expected for a system with global Z 2 symmetry. 

Methods. We outline the idea of fictitious vertices that we 
introduced in the current studyin the framework of the di¬ 
rected loop algorithm (DLA)^^“^^. We used the loop algo- 
rithm^^’'^^ in the basis diagonalizing erf’s for large Jxx supple¬ 
mented by cluster-type updates for elementary cubes in 3D. 
However, the loop algorithm leads to huge clusters in the QSL 
regime despite the very short correlation length, which re¬ 
duces the efficiency of the simulation. While such a difficulty 
is usually solved by the DLA, the conventional DLA encoun¬ 
ters an intrinsic ergodicity breakdown in the present model 
due to the off-diagonal multispin vertices. The fictitious ver¬ 
tex method can recover the broken ergodicity. This method 
is not necessarily a problem-specific idea, but a rather generic 
one"^^ to handle models with off-diagonal multispin vertices, 
such as spin-orbital"^^ or ring-exchange terms"^^""^^. 

We work on the cr^ basis, as this is convenient to measure 
fluxes. In the DLA^^“^^, the world-line configurations con- 
tributing to (n, T 2 ) = Tr exp{-f3^)af^ (T2) 

are sampled by moving one of the external lines (the “worm”). 



(b) 




FIG. 5. (Color online) (a) Worm-hopping at a A a vertex [illus¬ 
trated in (2-1-1)D] not allowed in the conventional DLA because 
{UAhihiul^s + Mthihihiu) = 0- (b) Example of a local 
configuration that cannot be generated in the conventional DLA. (c) 
By promoting the vertex into ?l fictitious vertex, the hopping of the 
type (a) can be included; the configurations like (b) can be realized if 
the subsequent worm move coincides with the directed dashed line. 


say at (i 1, Ti ), while fixing the other at (^2 , ^2 ) • At each space- 
time point, there is either the identity or one of the vertex op¬ 
erators, i.e.. As + 1 , Bp, or a faj 1 (which we call Aa, 
Xb, and Jxx vertices, respectively, and the constants are ab¬ 
sorbed into Jf), resulting from an expansion of exp(—eitf) to 
0 (e) with e <^T ^ . Every off-diagonal event in a world-line 
configuration corresponds to a vertex operator and its nonzero 
matrix element. The worm, locally updating the configuration 
as it moves, may be scattered at a vertex if the resulting state 
corresponds to a nonzero matrix element. However, an issue 
here is that this conventional scheme prohibits a hopping at 
Aa vertices [Fig. 5(a)]. Consequently, it is impossible to gen¬ 
erate off-diagonal processes that can arise from combinations 
of Aa and Jxx vertices [Fig. 5(b)]. This causes an intrinsic 
ergodicity issue in the conventional algorithm. 

Fictitious vertices allow such otherwise prohibited hopping 
by extending the configuration space. Similar ideas were dis¬ 
cussed in Refs. 48 and 49. For example, if the worm (at site 
1) reaches a A a vertex (at s), we promote this into a ficti¬ 
tious vertex {As + l)afaf, with a probability 1 — p, where 

(t^ 1) G 5 is chosen randomly (with a probability p, the 
worm follows the conventional scheme). With the aid of 
the attached a^ operators marking I and I', the worm hops 
from / to /' and leaves the vertex in ±r directions with equal 
probabilities [Fig. 5(c)]. This hopping process can be made 
rejection-free by setting the arbitrary weight for a fictitious 
vertex as Aa = (1 — p){zs — 1)~^Aa where = 4 ( 6 ) in 
2D (3D). After creating a fictitious vertex and the subsequent 
random-walk in the space-time^^, the worm may return to this 
fictitious vertex at a site yet to be marked by additional a ^. In 
2D, where As is a four-spin operator, we can bring this vertex 
back to normal after the worm hops to the last unmarked site. 
In 3D, after a similar process we still need to wait for the worm 
to return once again to recover a normal A a vertex. With such 
sequences we can generate off-diagonal multispin processes 
like the one in Fig. 5(b) or even more complicated ones. The 
observables are estimated when there is no fictitious vertex, 
namely, while we are sampling physical world-line configura¬ 
tions. 

Conclusions. We presented thermodynamic phase dia¬ 
grams of the 2D and 3D toric codes extended by a nearest- 
neighbor FM Ising coupling, which introduces quantum fluc¬ 
tuations in the magnetic flux excitations of the Z 2 QSLs, and 
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eventually stabilizes the FM phase. In 3D, we found a first- 
order transition between the Z 2 QSL and FM phases, with 
strong proximity effects to QSL in the phase competing region. 
The finite-T transition from paramagnet to QSL is in the in¬ 
verted Ising universality class, as in the case without the Ising 
coupling. These characteristics will provide a clue for exper¬ 
imental studies on the 3D Z 2 QSL candidates. In 2D, the Z 2 
QSL state is unstable at finite T and forms a QCP between the 
competing FM phase. We found the quantum critical behavior 
in the (2-1-1 )D Ising universality class with an unconventional 
exponent for the magnetic susceptibility. The distinct behav¬ 
iors between 2D and 3D are essentially due to whether flux 
excitations are deconfined particles (in 2D) or confined loops 


(in 3D) at low T. 
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